Trabajo Práctico 1: Regresión lineal
Mauro Bringas
2023-10-19
Importamos librerias
library(recipes)## Loading required package: dplyr##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union##
## Attaching package: 'recipes'## The following object is masked from 'package:stats':
##
## steplibrary(corrplot)## corrplot 0.92 loadedlibrary(corrr)
library(dplyr)
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats 1.0.0 ✔ readr 2.1.4
## ✔ ggplot2 3.4.4 ✔ stringr 1.5.0
## ✔ lubridate 1.9.2 ✔ tibble 3.2.1
## ✔ purrr 1.0.2 ✔ tidyr 1.3.0## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ stringr::fixed() masks recipes::fixed()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(rsample)
library(GGally)## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2library(robust)## Loading required package: fit.modelslibrary(magrittr)##
## Attaching package: 'magrittr'
##
## The following object is masked from 'package:purrr':
##
## set_names
##
## The following object is masked from 'package:tidyr':
##
## extractrequire(generics)## Loading required package: generics
##
## Attaching package: 'generics'
##
## The following object is masked from 'package:lubridate':
##
## as.difftime
##
## The following object is masked from 'package:dplyr':
##
## explain
##
## The following objects are masked from 'package:base':
##
## as.difftime, as.factor, as.ordered, intersect, is.element, setdiff,
## setequal, unionrequire(yardstick)## Loading required package: yardstick
##
## Attaching package: 'yardstick'
##
## The following object is masked from 'package:generics':
##
## accuracy
##
## The following object is masked from 'package:readr':
##
## spec1) Análisis de estructura y correlación
Carga de datos
eph_train0 <- read.csv("eph_train_2022.csv")Leer el dataset
Leer el archivo “eph_train_2022.csv”. ¿Qué puede mencionar sobre su estructura y variables?
#Vemos un poco de los datos
glimpse(eph_train0)## Rows: 11,625
## Columns: 18
## $ codusu <chr> "TQRMNOQTQHJOKRCDEHPJB00784929", "TQRMNOPRTHKMLS…
## $ ano4 <int> 2022, 2022, 2022, 2022, 2022, 2022, 2022, 2022, …
## $ trimestre <int> 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, …
## $ region <chr> "Noroeste", "Pampeana", "Noroeste", "Noroeste", …
## $ aglomerado <int> 29, 13, 19, 19, 10, 3, 38, 2, 36, 26, 8, 19, 10,…
## $ fecha_nacimiento <chr> "28/08/1961", "01/01/1900", "27/01/1982", "11/11…
## $ edad <int> 60, 60, 40, 21, 62, 20, 64, 39, 21, 39, 44, 51, …
## $ asistencia_educacion <chr> "asistio", "asistio", "asistio", "asistio", "asi…
## $ nivel_ed <chr> "Primaria\nIncompleta", "Superior\nUniversitaria…
## $ tipo_establecimiento <chr> "Privada", "Privada", "Estatal", "Estatal", "Pri…
## $ codigo_actividad <int> 3100, 4803, 8600, 8401, 6100, 4803, 4000, 8600, …
## $ sexo <chr> "Varon", "Varon", "Mujer", "Varon", "Varon", "Va…
## $ categoria_ocupacion <chr> "Cuenta propia", "Cuenta propia", "Obrero o empl…
## $ cat_cantidad_empleos <chr> "unico", "unico", "unico", "unico", "unico", "un…
## $ alfabetismo <chr> "Sabe leer y escribir", "Sabe leer y escribir", …
## $ educacion <int> 6, 15, 13, 12, 15, 13, 4, 19, 13, 8, 19, 16, 16,…
## $ experiencia_potencial <int> 49, 40, 22, 4, 42, 2, 55, 15, 3, 26, 20, 30, 22,…
## $ salario_horario <dbl> 214.28571, 83.33333, 541.66667, 500.00000, 619.3…El dataset tiene 18 variables. Un código de identificacion, el año y trimestre al cual corresponde el registro son atributos identificatorios de este tipo de datasets verticales (donde cada registro es un punto temporal definido para un dado individuo).
La region, variable de tipo caracter, y el aglomerado, numérica, hacen referencia a la ubicacion del caso registrado.
La fecha de nacimiento, edad, asistencia al sistema educativo, el último nivel educativo alcanzado, el alfabetismo y el sexo son caracteristicas demográficas.
En tercera instancia tenemos características de su actividad laboral, como ser el código de actividad, la categoria de la ocupación, la experiencia laboral y el salario horario.
En este trabajo el salario horario será nuestra variable objetivo. Nos sacaremos de encima las variables identificatorias
eph_train <- eph_train0 %>% select(-c(codigo_actividad,codusu,aglomerado,ano4,trimestre))Evaluamos si hay valores faltantes
#Seleccionamos variables numericas y passamos a dataset en formato largo. Vemos si hay valores perdidos
eph_train_vert <- eph_train %>% select(where(is.numeric),sexo) %>% pivot_longer(.,col=-c(sexo),names_to = "variables",values_to = "valores") %>% select(-c(sexo))
eph_train_vert %>% gather(.,
key = "variables",
value = "valores") %>%
group_by(variables) %>%
summarise(valores_unicos = n_distinct(valores),
porcentaje_faltantes = sum(is.na(valores))/nrow(eph_train_vert)*100) %>%
arrange(desc(porcentaje_faltantes), valores_unicos) ## # A tibble: 4 × 3
## variables valores_unicos porcentaje_faltantes
## <chr> <int> <dbl>
## 1 educacion 24 0
## 2 experiencia_potencial 73 0
## 3 edad 74 0
## 4 salario_horario 1623 0En ninguna de las variables numéricas de interés hay valores faltantes.
eph_train_vert <- eph_train %>% select(-where(is.numeric),sexo) %>% pivot_longer(.,col=-c(sexo),names_to = "variables",values_to = "valores") %>% select(-c(sexo))
eph_train_vert %>% gather(.,
key = "variables",
value = "valores") %>% # agrupamos por las variables del set
group_by(variables) %>%
summarise(valores_unicos = n_distinct(valores),
porcentaje_faltantes = sum(is.na(valores))/nrow(eph_train_vert)*100) %>%
arrange(desc(porcentaje_faltantes), valores_unicos) ## # A tibble: 8 × 3
## variables valores_unicos porcentaje_faltantes
## <chr> <int> <dbl>
## 1 alfabetismo 2 0
## 2 asistencia_educacion 2 0
## 3 cat_cantidad_empleos 2 0
## 4 tipo_establecimiento 3 0
## 5 categoria_ocupacion 4 0
## 6 nivel_ed 6 0
## 7 region 6 0
## 8 fecha_nacimiento 7428 0Tampoco esto ocurre en las variables categóricas.
¿Cómo es la correlación entre las variables numéricas? Utilice y analice en detalle algún gráfico que sirva para sacar conclusiones sobre la asociación de variables realizando apertura por sexo. En particular, ¿Cómo es la correlación entre la variable a explicar (salario_horario) y el resto de las variables numéricas?
eph_train %>%
select(where(is.numeric),sexo) %>% # desestimamos algunas variables
mutate(sexo = factor(sexo)) %>%
ggpairs(., mapping = aes(colour = sexo), title = "Matriz de correlaciones",
upper = list(continuous = wrap("cor", size = 3, hjust=0.5)), legend = 25) +
theme_bw() +
theme(axis.text.x = element_text(angle=45, vjust=0.5), legend.position = "bottom")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Podemos ver que individualmente tanto la edad como la experiencia potencial y los años de educacion correlacionan positivamente con el salario horario con coeficientes 0.16, 0.06 y 0.35, respectivamente.
2) Modelo en funcion de la experiencia laboral y de la experiencia laboral al cuadrado.
Se va a comenzar con dos modelos lineales que utilicen la información de la experiencia potencial. Primero, ajustar un modelo de regresión para explicar el salario por hora usando únicamente la experiencia potencial como covariable.
\(E(SalarioHorario) = β0 + β1*ExperienciaPotencial\)
Luego, ajustar otro modelo en donde las únicas covariables sean la experiencia potencial y el cuadrado de la experiencia potencial.
\(E(SalarioHorario) = β0 + β1*ExperienciaP otencial + β2*ExperienciaPotencial^2\)
Responder las siguientes preguntas en base a ambos modelos: ¿Cuál es el impacto de un año adicional de experiencia potencial en el salario horario esperado para cada uno de estos modelos? ¿Cuál es el efecto sobre el salario horario esperado de un año más de experiencia laboral para una persona con 6 años de experiencia laboral? ¿Y para una persona con 35 años de experiencia laboral?
mls_exp <- eph_train %>% lm(salario_horario~experiencia_potencial,.)
summary(mls_exp)##
## Call:
## lm(formula = salario_horario ~ experiencia_potencial, data = .)
##
## Residuals:
## Min 1Q Median 3Q Max
## -581.4 -250.9 -98.2 111.2 4448.9
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 480.1808 7.4701 64.28 < 2e-16 ***
## experiencia_potencial 1.8664 0.2841 6.57 5.23e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 410.9 on 11623 degrees of freedom
## Multiple R-squared: 0.003701, Adjusted R-squared: 0.003615
## F-statistic: 43.17 on 1 and 11623 DF, p-value: 5.23e-11En este modelo la ordenada al origen indica que, en el caso de una persona sin experiencia potencial, el salario horario medio esperado es 480.18 pesos. Por otro lado, la pendiente sugiere que por cada año de experiencia laboral incrementada, el salario horario medio se incrementa en 1.87 pesos.
#Obtenemos los coeficientes y sus intervalos de confianza para graficar
tidy_model <- broom::tidy(mls_exp,conf.int = TRUE)
ggplot(tidy_model, aes(estimate, term, xmin = conf.low, xmax = conf.high, height = 0)) +
geom_point(color = "forestgreen",size=2) +
geom_vline(xintercept = 0, lty = 4, color = "black") +
geom_errorbarh(color = "forestgreen", size=1) +
theme_bw() +
labs(y = "Coeficientes β", x = "Estimación")## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
Podemos ver que tanto la ordenada al origen como la pendiente tienen p-valores muy pequeños, asegurando que el modelo realizado es significativamente distinto de modelar este problema con una constante (y esa constante es distinta de cero). De forma análoga podemos decir que, por ejemplo, a nivel de confianza 0.95 los intervalos alrededor de los coeficientes estimados no incluyen al cero.
Al ser una regresion lineal simple tenemos la ventaja de que podemos graficar la variable a predecir en funcion de la covariable
#Tomamos pendiente y ordenada al origen
intercepto = mls_exp$coefficients[1]
pendiente = mls_exp$coefficients[2]
#Graficamos los puntos y el ajuste con una recta
eph_train %>% ggplot(., aes(x = experiencia_potencial, y = salario_horario)) +
geom_point() +
geom_abline(intercept = intercepto, slope = pendiente, color="forestgreen") +
labs(title="Modelo Lineal Simple", x="Experiencia Potencial", y="Salario Horario") 
Ahora incorporamos el cuadrado de la experiencia potencial a nuestro modelo lineal simple.
mls_exp_exp2 <- eph_train %>% mutate(exp2=experiencia_potencial**2) %>% lm(salario_horario~experiencia_potencial+exp2,.)
summary(mls_exp_exp2)##
## Call:
## lm(formula = salario_horario ~ experiencia_potencial + exp2,
## data = .)
##
## Residuals:
## Min 1Q Median 3Q Max
## -539.9 -248.2 -100.3 113.1 4443.3
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 387.72725 11.12290 34.86 <2e-16 ***
## experiencia_potencial 12.16103 0.96395 12.62 <2e-16 ***
## exp2 -0.20301 0.01817 -11.17 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 408.7 on 11622 degrees of freedom
## Multiple R-squared: 0.01428, Adjusted R-squared: 0.01411
## F-statistic: 84.2 on 2 and 11622 DF, p-value: < 2.2e-16El resultado del modelo lineal otorga tres parametros cuya significacion se corrobora con un p valor menor a 10^-15.
#Obtenemos los coeficientes y sus intervalos de confianza para graficar
tidy_model <- broom::tidy(mls_exp_exp2,conf.int = TRUE)
ggplot(tidy_model, aes(estimate, term, xmin = conf.low, xmax = conf.high, height = 0)) +
geom_point(color = "forestgreen",size=2) +
geom_vline(xintercept = 0, lty = 4, color = "black") +
geom_errorbarh(color = "forestgreen", size=1) +
theme_bw() +
labs(y = "Coeficientes β", x = "Estimación")
La ordenada al origen indica que el salario por hora esperado para los individuos sin experiencia laboral es de 387.72. Por otro lado, el término que acompaña a la experiencia potencial indica que por cada año dicional de experiencia potencial el salario por hora esperado se incrementa en 12.16 pesos. El término cuadrático (en experiencia potencial) en cambio tiene signo negativo, lo cual indica que existe un punto de inflexión para el comportamiento de la experiencia potencial. Este término es el análogo a la aceleración si hacemos una comparación entre tiempo-experiencia y salario por hora - altura en un tiro vertical.
Para determinar cuál es el punto de inflexion pasamos esa expresión a la forma canónica de los polinomios de grado 2.
\(SalarioHorario = -0.20301 * (ExperienciaPotencial-29.952)^2 + 69.85\)
Esto muestra que el crecimiento de la esperanza de salario horario con la experiencia será cada vez más grande, hasta que en los 29.952 años de experiencia esa tasa de variación comenzará a decrecer.
Tambien podriamos pensar que en este modelo lineal la constante que acompaña a la experiencia potencial tiene una dependencia con esa misma variable
\(E(SalarioHorario) = β_0 + β_X*ExperienciaPotencial\)
\(β_X = β_1 + β_2*ExperienciaPotencial\)
Reiterando la pregunta del enunciado, ¿Cuál es el efecto sobre el salario horario esperado de un año más de experiencia laboral para una persona con 6 años de experiencia laboral? ¿Y para una persona con 35 años de experiencia laboral?
Miramos entonces la variación respecto del salario potencial (o sea, la derivada primera de la esperanza en salario horario respecto de la experiencia potencial)
\(\frac{d_{salario}}{d_{experienciapotencial}}=12.16103-2*0.20301*experienciaPotencial\)
En una persona con 6 años, aumentar un año aumenta 9.72 pesos/hora el valor esperado del salario horario y en una con 35 años, luego del punto de inflexión, la esperanza del salario horario tiene una variacion de -2.05 pesos/hora.
#Graficamos los puntos y el ajuste polinomico de orden 2 (o sea, el modelo lineal con features experiencia y experiencia al cuadrado)
eph_train %>% ggplot(., aes(x = experiencia_potencial, y = salario_horario)) +
geom_point() + stat_smooth(method='lm', formula = y~poly(x,2),se=FALSE)+
labs(title="Modelo Lineal Simple con Experiencia al cuadrado", x="Experiencia Potencial", y="Salario Horario") 
¿Cuál es la interpretación de las variables incluidas en el modelo? ¿Sus coeficientes son significativos? ¿El modelo resulta significativo para explicar el salario? ¿Qué porcentaje de la variabilidad explica el modelo?
En el modelo simple que solo incluye la experiencia laboral linealmente, la ordenada al origen es el valor esperado de salario horario para una persona sin experiencia potencial. El coeficiente que acompaña a la experiencia potencial representa la variacion en el valor medio esperado de salario horario.
En el modelo que incluye la experiencia laboral al cuadrado, la interpretación de ambos coeficientes (que acompañan a experiencia potencial y experiencia potencial al cuadrado) debe realizarse al mismo tiempo, dado que estas covariables se mueven de forma totalmente correlacionada. De esta forma, podemos tomar en consideracion que la variacion en experiencia potencial no es igual para cualquier valor de experiencia potencial.
3) Modelo lineal múltiple
Se plantea un primer modelo múltiple a partir de la ecuación de Mincer: \(E(SalarioHorario) = β_0 + β_1 AñosEducacion + β_2 ExperienciaPotencial +\) \(β_2 ExperienciaP otencial^2 + β_3 Sexo + β_4 Sexo·AñosEducacion\)
Ajustar el modelo planteado y responder las siguientes preguntas: ¿Cuál es la interpretación de las variables incluidas en el modelo? ¿Sus coeficientes son significativos? ¿El modelo resulta significativo para explicar el salario? ¿Qué porcentaje de la variabilidad explica el modelo?. Analizar en profundidad el cumplimiento de los supuestos del modelo lineal para este modelo.
En primera instancia ajustamos el modelo
mlm <- eph_train %>% mutate(exp2=experiencia_potencial**2) %>% lm(salario_horario~educacion+experiencia_potencial+exp2+sexo+sexo*educacion,.)
summary(mlm)##
## Call:
## lm(formula = salario_horario ~ educacion + experiencia_potencial +
## exp2 + sexo + sexo * educacion, data = .)
##
## Residuals:
## Min 1Q Median 3Q Max
## -876.4 -211.2 -66.6 114.8 4092.6
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -319.8261 24.3530 -13.133 < 2e-16 ***
## educacion 48.2334 1.5115 31.911 < 2e-16 ***
## experiencia_potencial 11.1425 0.8889 12.535 < 2e-16 ***
## exp2 -0.1019 0.0169 -6.031 1.68e-09 ***
## sexoVaron 73.2669 26.9747 2.716 0.00661 **
## educacion:sexoVaron -2.0333 1.9573 -1.039 0.29891
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 376.6 on 11619 degrees of freedom
## Multiple R-squared: 0.1634, Adjusted R-squared: 0.163
## F-statistic: 453.7 on 5 and 11619 DF, p-value: < 2.2e-16En este caso podemos ver que el p valor del test F para el modelo global da que es significativo, lo cual dice que realizar este modelo difiere significativamente de modelar el salario horario con una constante.
Interpretación:
La introduccion de una variable categorica hace que podamos pensar modelos “distintos” para cada grupo. En este caso, al estar indicado sexoVaron en los coeficientes, notamos que la categoria de referencia es Mujer. Por esto, el modelo de las Mujeres tendra valor sexoVaron=0, y el modelo de varones sexoVaron=1.
La ordenada al origen es el salario horario esperado para un individuo sin experiencia laboral, de sexo femenino y cero años de educacion informados.
El coeficiente asociado a sexoVaron indica cuál es la diferencia en salario horario esperado para un individuo sin experiencia potencial ni educacion, considerando que su sexo es Varon.
El coeficiente que acompaña a educacion representa el incremento en salario horario esperado por cada año de educacion incrementado, para un individuo de sexo femenino y a igualdad de condiciones de experiencia potencial.
El coeficiente que acompaña al termino de interaccion educacion:sexoVaron representa cuánto más se modifica el sarlario horario esperado por cada año de educacion incrementado por encima de lo dicho para el termino que afecta a individuos de sexo femenino, dejando constante la experiencia potencial.
El coeficiente que acompaña a experiencia y experiencia al cuadrado muestra cómo una variacion en experiencia potencial repercute sobre el salario horario esperado, dejando constante el resto de las covariables.
No todos los coeficientes descritos son significativos (no existe evidencia suficiente para afirmar que son distintos de cero, a un nivel de confianza de 0.95), en particular educacion*sexoVaron tiene un pvalor de casi 0.3. Esto quiere decir que no existe suficiente evidencia para afirmar que la influencia de incrementar los años de educacion en varones, por sobre las mujeres, sea distinto de cero para cualquier nivel de confianza mayor que 70%.
#Obtenemos los coeficientes y sus intervalos de confianza para graficar
tidy_model <- broom::tidy(mlm,conf.int = TRUE)
ggplot(tidy_model, aes(estimate, term, xmin = conf.low, xmax = conf.high, height = 0)) +
geom_point(color = "forestgreen",size=2) +
geom_vline(xintercept = 0, lty = 4, color = "black") +
geom_errorbarh(color = "forestgreen", size=1) +
theme_bw() +
labs(y = "Coeficientes β", x = "Estimación")
Pasando en limpio, el modelo de las mujeres es
\(E(SalarioHorario) = β_0 + β_1 AñosEducacion + β_2 ExperienciaPotencial + β_2 ExperienciaPotencial^2\)
y el de los varones,
\(E(SalarioHorario) = (β_0 + β_3) + (β_1 + β_4) AñosEducacion + β_2 ExperienciaPotencial + β_2 ExperienciaP otencial^2\)
En esta instancia vale la pena explicar qué implica que sea tan grande el p-valor del coeficiente β_4 (educacion para el sexo Varon): No es significativa la diferencia en la influencia de los años de educacion de los varones por sobre lo que se estima para las mujeres en el valor esperado del salario horario, a todas las otras variables fijas (en otras palabras, no existe suficiente evidencia para afirmar que β_4 sea distinto de cero).
Significatividad global del modelo
Podemos observar en la salida del ajuste del modelo lineal que el valor observado para el estadístico F se corresponde con un p-valor muy pequeño, menor a \(2.10^{-16}\). Esto quiere decir que globalmente el modelo es útil para modelar el salario horario esperado en función de las covariables utilizadas, más allá de que individualmente algunos coeficientes hayan resultado no significativos a niveles de confianza relativamente bajos.
Explicacion de la varianza
La fracción de varianza explicada en los modelos lineales se obtiene del coeficiente \(R^2\), Al introducir más de una covariable, se realiza un ajuste de este coeficiente para tomar en consideración que en general al aumentar indefinidamente las covariables, podría incrementarse la fracción de varianza explicada pero de manera estocástica y no explicativa. El coeficiente \(R^2_{adj}\) en este caso es de 0.163, lo cual indica aproximadamente un 16% de la varianza total explicada. Cabe resaltar que el \(R^2\) sin ajustar devuelve una varianza explicada similar, lo cual habla sobre un aporte apreciable de las covariables utilizadas para la explicación de la varianza.
Cumplimiento de supuestos
Los supuestos que debemos verificar son
-Homocedasticidad de los residuos
-Normalidad de los residuos
-Independencia de los residuos (lo cual queda determinado por el diseño experimental)
Para esto, utilizamos la función plot sobre el modelo lineal generado y así obtener un grafico de residuos versus valores predichos, un Q-Q Plot, un gráfico de ubicación extendida y uno de residuos versus potencial de apalancamiento.
#Obtenemos los graficos para hacer el diagnóstico.
plot(mlm)
Del grafico de residuos versus valores predichos podemos ver que parece existir cierta estructura en los datos dado que hay una curvatura. A bajos y altos valores los residuos son positivos, mientras que en valores intermedios son negativos.
Del Q-Q Plot de residuos estandarizados podemos ver que la distribucion se aleja bastante de la comparacion cuantil-a-cuantil con la distribucion N(0,1). Esto es Particularmente notable a valores de cuantiles teóricos altos y muestra una distribución asimetrica y tirada hacia la izquierda. Esto viola el supuesto de normalidad de los residuos.
Del grafico de residuos versus potencial apalancamiento podemos afirmar que no existen puntos potencialmente influyentes, dado que 0.5>0.012>hii observados.
El diagnóstico para este modelo es: modelo no cumple con los supuestos del modelo lineal. Parecen existir dos problemas: violación del supuesto de linealidad de la esperanza condicional (residuos con estructura), falta de normalidad en los residuos(del Q-Q Plot).
4) Modelo de Mincer “enriquecido”
Ahora, se procede a modelar según una especificación del modelo de Mincer con ciertas variables adicionales
\(E[ln(SalarioHorario)] = β0 + β1AñosEducacion + β2ExperienciaP otencial + β3ExperienciaP otencial2+ β4Sexo + β5Sexo · AñosEducacion\)
• ¿Cuál es la interpretación del coeficiente asociado a la variable de años de educación? ¿Se observan cambios en la significatividad individual de los coeficientes respecto al modelo anterior?
• ¿Qué porcentaje de la variabilidad del salario horario explica el modelo? ¿Cómo se compara con la variabilidad explicada por el modelo anterior?
Nota: tenga en cuenta que la variable predicha es el logaritmo del salario horario y se pide el porcentaje de variabilidad explicada del salario horario. Además, como los dos modelos tienen la misma cantidad de covariables es posible compararlos mediante el el R-cuadrado simple.
• Analizar en profundidad el cumplimiento de los supuestos del modelo lineal para este modelo y comparar con el análisis del modelo anterior.
lmMincerEnriquecido <- eph_train %>% mutate(exp2=experiencia_potencial**2) %>% lm(log(salario_horario) ~ educacion+experiencia_potencial+exp2 + sexo +
sexo*educacion,.)
summary(lmMincerEnriquecido)##
## Call:
## lm(formula = log(salario_horario) ~ educacion + experiencia_potencial +
## exp2 + sexo + sexo * educacion, data = .)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.5868 -0.3629 0.0399 0.3950 2.8591
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.441e+00 4.083e-02 108.760 < 2e-16 ***
## educacion 8.969e-02 2.534e-03 35.391 < 2e-16 ***
## experiencia_potencial 2.366e-02 1.490e-03 15.879 < 2e-16 ***
## exp2 -2.774e-04 2.834e-05 -9.790 < 2e-16 ***
## sexoVaron 2.480e-01 4.522e-02 5.483 4.27e-08 ***
## educacion:sexoVaron -1.141e-02 3.282e-03 -3.478 0.000507 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6314 on 11619 degrees of freedom
## Multiple R-squared: 0.183, Adjusted R-squared: 0.1827
## F-statistic: 520.6 on 5 and 11619 DF, p-value: < 2.2e-16La variable a modelar en este caso es el logaritmo del salario horario, y todo esto constituye un modelo semi elastico. En esta oportunidad, todos los coeficientes asociados a las distintas covariables resultan significativos incluso a niveles de confianza de 0.999.
#Obtenemos los coeficientes y sus intervalos de confianza para graficar
tidy_model <- broom::tidy(lmMincerEnriquecido,conf.int = TRUE)
ggplot(tidy_model, aes(estimate, term, xmin = conf.low, xmax = conf.high, height = 0)) +
geom_point(color = "forestgreen",size=2) +
geom_vline(xintercept = 0, lty = 4, color = "black") +
geom_errorbarh(color = "forestgreen", size=1) +
theme_bw() +
labs(y = "Coeficientes β", x = "Estimación")
Ningun intervalo de confianza abarca al valor cero.
Al tener un modelo semielastico (Modelo Log-Nivel) -la variable a explicar en forma logarítmica y las covariables en forma lineal- la interpretación del coeficiente se realiza de la siguiente manera: en particular para el caso de años de educacion, por cada año de incremento de la variable educacion, se tiene un cambio esperado de 8,97% el salario horario.
Porcentaje de variabilidad explicado
Una comparación de \(R^2\) será posible sólo si la variable dependiente es la misma para ambos modelos. En este caso, el modelo lineal multiple y el semielastico utilizan las mismas covariables, con lo cual es posible comparar el porcentaje de explicación de la varianza.
El \(R^2\) simple y ajustado en ambos modelos no difiere mucho en cada caso, con lo cual la comparacion entre modelos la haremos con el \(R^2\) simple. Para el modelo Nivel-Nivel es de 0.1634 y en el modelo analizado en este item del TP es de 0.1763. Esto quiere decir que ambos modelos explican entre el 16 y 17% de la varianza total, observando un incremento para el modelo Log-Nivel.
Cumplimiento de supuestos
Ahora analizaremos el cumplimiento de los supuestos del modelo lineal para este modelo y compararemos con el análisis del modelo lineal múltiple Nivel-Nivel.
#Obtenemos los graficos para realizar el diagnóstico
plot(lmMincerEnriquecido)
En este caso, los residuos en funcion de los valores predichos siguen teniendo estructura (siendo mayormente positivos a altos y bajos valores de valores predichos).
Del grafico del Q-Q Plot podemos observar que tampoco se cumple el supuesto de normalidad, pero en este caso las desviaciones en cuantiles bajos y altos son distintas. En este caso, la distribucion de residuos tiene las colas pesadas comparadas con la distribución normal.
En el grafico de potencial apalancamiento tampoco observamos valores que pudieran llegar a generar problemas por tener alto potencial de apalancamiento y ubicarse en una región atípica.
Ambos modelos tienen los mismos problemas respecto de los residuos calculados. No cumplen la homocedasticidad ni normalidad.
5) Modelos propios y evaluación
Realizar 2 modelos lineales múltiples adicionales y explicar la lógica detrás de los mismos (se valorará la creación y/o inclusión de variables nuevas).
Nota: No se pueden utilizar métodos de selección automática de variables dado que buscamos que analicen otras variables y realicen feature engineering.
Evaluar y comparar la performance del modelo lineal multiple, el modelo de mincer y los modelos desarrollados en este punto en el dataset de entrenamiento y evaluación (usar dataset “eph_test_2022.csv”).
La evaluación de performance consiste en comparar la performance en términos del RMSE y MAE sobre el set de entrenamiento y el set de evaluación.
¿Cuál es el mejor modelo para el objetivo de predecir el salario horario? ¿Por qué?
El primer modelo propuesto involucra la creacion de una nueva variable categorica, que junta varias regiones en distintas zonas, y será incorporado al modelo únicamente afectando al término independiente y no generando términos de interacción con otras variables.
#Creamos la nueva variable categorica
lmPropio1 <- eph_train %>%
mutate(zona = case_when(region == "Capital Federal" ~ "AMBA",
region == "Gran Buenos Aires" ~ "AMBA",
region == "Noreste" ~ "Norte",
region == "Noroeste" ~ "Norte",
region == "Patagonia" ~ "Sur",
region == "Cuyo" ~ "Centro",
region == "Pampeana" ~ "Centro"),
exp2=experiencia_potencial^2) %>%
lm(salario_horario ~ educacion+experiencia_potencial+exp2 +zona,.)
summary(lmPropio1)##
## Call:
## lm(formula = salario_horario ~ educacion + experiencia_potencial +
## exp2 + zona, data = .)
##
## Residuals:
## Min 1Q Median 3Q Max
## -853.3 -200.9 -61.6 106.6 4116.6
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -164.13782 19.92902 -8.236 < 2e-16 ***
## educacion 45.19502 1.00056 45.170 < 2e-16 ***
## experiencia_potencial 10.84489 0.86868 12.484 < 2e-16 ***
## exp2 -0.10166 0.01653 -6.150 8.02e-10 ***
## zonaCentro -76.43465 11.05371 -6.915 4.93e-12 ***
## zonaNorte -175.71263 11.18410 -15.711 < 2e-16 ***
## zonaSur 71.92236 13.53976 5.312 1.10e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 367.9 on 11618 degrees of freedom
## Multiple R-squared: 0.2017, Adjusted R-squared: 0.2013
## F-statistic: 489.3 on 6 and 11618 DF, p-value: < 2.2e-16Este modelo incorpora un agrupamiento de zonas geograficas (Centro, Norte, Sur, AMBA) y toma al AMBA como referencia. En este caso se introdujo esta covariable discreta llamada “Zona” sin interaccion con otros términos, proponiendo que la influencia de la educacion y la experiencia laboral en todas las regiones será similar en todas las zonas. Bajo este modelo, solo se modificará el salario horario a educacion y experiencia laboral constante, al moverse de región a región.
#Obtenemos los coeficientes y sus intervalos de confianza para graficar
tidy_model <- broom::tidy(lmPropio1,conf.int = TRUE)
ggplot(tidy_model, aes(estimate, term, xmin = conf.low, xmax = conf.high, height = 0)) +
geom_point(color = "forestgreen",size=2) +
geom_vline(xintercept = 0, lty = 4, color = "black") +
geom_errorbarh(color = "forestgreen", size=1) +
theme_bw() +
labs(y = "Coeficientes β", x = "Estimación")
En este caso, todos los coeficientes son significativamente distintos de cero y pueden interpretarse de la siguiente manera:
La ordenada al origen es el valor esperado de salario minimo para un individuo en el AMBA con 0 años de educacion y 0 años de experiencia potencial.
El coeficiente de cada zona (Centro, Norte, Sur) indica cuanto difiere el salario horario esperado, a educacion y experiencia constante, al tomar un individuo y pensarlo en otra zona que no sea el AMBA. Por ejemplo zonaNorte dice que el salario horario esperado a igual educacion y experiencia potencial será 175.71 pesos/hora menor que el de el AMBA.
La interpretación de educacion y experiencia potencial es análoga a los modelos analizados anteriormente.
Para un segundo modelo se propone incorporar una variable continua que muestre el porcentaje de vida en el que un individuo ha acumulado experiencia, variando entre 0 (no tiene experiencia) y 1 (trabajó todos los años de su vida, caso hipotético).
#Creamos la nueva variable numerica
lmPropio2 <- eph_train %>%
mutate(vida_trabajada = experiencia_potencial/edad,
exp2=experiencia_potencial^2) %>%
lm(salario_horario ~ educacion+experiencia_potencial+exp2+
vida_trabajada,.)
summary(lmPropio2)##
## Call:
## lm(formula = salario_horario ~ educacion + experiencia_potencial +
## exp2 + vida_trabajada, data = .)
##
## Residuals:
## Min 1Q Median 3Q Max
## -882.1 -212.5 -65.5 115.9 4127.3
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -546.7066 44.1945 -12.370 < 2e-16 ***
## educacion 59.0425 2.1226 27.817 < 2e-16 ***
## experiencia_potencial -16.0352 3.9598 -4.049 5.17e-05 ***
## exp2 0.1847 0.0440 4.199 2.71e-05 ***
## vida_trabajada 1044.4557 148.9929 7.010 2.51e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 376.5 on 11620 degrees of freedom
## Multiple R-squared: 0.1638, Adjusted R-squared: 0.1635
## F-statistic: 569 on 4 and 11620 DF, p-value: < 2.2e-16En este caso nuevamente todos los coeficientes son significativos, y el agregado distintivo de la variable creada vida_Trabajada indica cómo se modifica el salario horario esperado a igualdad de educacion y experiencia potencial, cuando cambia la fraccion de vida en la cual se acumuló experiencia potencial.
En este caso, este coeficiente afirma que a igualdad de educación y experiencia potencial, una persona con cociente entre experiencia y edad muy alto (tendiendo a uno) tendra un salario horario esperado 1044.45 pesos/hora más alto que una persona donde este cociente sea muy pequeño (tendiente a cero).
#Obtenemos los coeficientes y sus intervalos de confianza para graficar
tidy_model <- broom::tidy(lmPropio2,conf.int = TRUE)
ggplot(tidy_model, aes(estimate, term, xmin = conf.low, xmax = conf.high, height = 0)) +
geom_point(color = "forestgreen",size=2) +
geom_vline(xintercept = 0, lty = 4, color = "black") +
geom_errorbarh(color = "forestgreen", size=1) +
theme_bw() +
labs(y = "Coeficientes β", x = "Estimación")
En este caso también todos los coeficientes resultan significativos, incluso a nivel de confianza 0.999 .
Evaluacion de modelos
Haremos un análisis de cumplimiento de supuestos en ambos modelos llamamdos lmPropio1 y lmPropio2.
Comenzamos con el modelo propio 1:
#Obtenemos los graficos para realizar el diagnóstico
plot(lmPropio1)
Podemos observar en el gráfico de residuos versus valores predichos que la varianza de los residuos parece incrementarse a medida que crece el valor predicho. No se cumple la homocedasticidad.
En segundo lugar, el Q-Q Plot muestra que los cuantiles de la distribucion observada no se corresponden con los de la normal. En este caso, la distribución es asimétrica a izquierda. No se cumple el supuesto de normalidad de los residuos.
Continuamos con el modelo propio 2:
#Obtenemos los graficos para realizar el diagnóstico
plot(lmPropio2)
Podemos observar en el gráfico de residuos versus valores predichos que la dispersion de los residuos crece a medida que aumenta el valor predicho. No se cumple la homocedasticidad.
En segundo lugar, el Q-Q Plot muestra que los cuantiles de la distribucion observada no se corresponden con los de la normal. Nuevamente, la distribución de resiudos es asimétrica a izquierda. No se cumple el supuesto de normalidad de los residuos.
Comparación de la performance del modelo lineal multiple, el modelo de mincer y los modelos desarrollados en este punto
Las dos métricas utilizadas para evaluar la performance son RMSE y MAE. Las calculamos para los cuatro modelos en entrenamiento y validación.
eph_test <- read.csv("eph_test_2022.csv")%>% select(-c(codigo_actividad,codusu,aglomerado,ano4,trimestre))
models <- list(mlm=mlm,lmMincerEnriquecido = lmMincerEnriquecido, lmPropio1= lmPropio1, lmPropio2 = lmPropio2)En primera instancia evaluamos las métricas en entrenamiento:
eph_train <- eph_train%>% mutate(exp2=experiencia_potencial**2,
zona = case_when(region == "Capital Federal" ~ "AMBA",
region == "Gran Buenos Aires" ~ "AMBA",
region == "Noreste" ~ "Norte",
region == "Noroeste" ~ "Norte",
region == "Patagonia" ~ "Sur",
region == "Cuyo" ~ "Centro",
region == "Pampeana" ~ "Centro"),
vida_trabajada = experiencia_potencial/edad
)
# Obtenemos las predicciones de ambos modelos
lista_predicciones_train = map(.x = models, .f = broom::augment)
metricas_mlm = lista_predicciones_train$mlm %>%
metrics(truth=salario_horario, estimate=.fitted) %>%
mutate(.estimate=round(.estimate, 4))
metricas_mlm$modelo <- "mlm"
metricas_lmPropio1 = lista_predicciones_train$lmPropio1%>%
metrics(truth=salario_horario, estimate=.fitted) %>%
mutate(.estimate=round(.estimate, 4))
metricas_lmPropio1$modelo <- "lmPropio1"
metricas_lmPropio2 = lista_predicciones_train$lmPropio2 %>%
metrics(truth=salario_horario, estimate=.fitted) %>%
mutate(.estimate=round(.estimate, 4))
metricas_lmPropio2$modelo <- "lmPropio2"
metricas_semilog = lista_predicciones_train$lmMincerEnriquecido %>%
mutate(fitted_antilog= exp(.fitted),salario_horario=exp(`log(salario_horario)`)) %>%
metrics(truth=salario_horario, estimate=fitted_antilog) %>%
mutate(.estimate=round(.estimate, 4))
metricas_semilog$modelo <- "lmMincerEnriquecido"
metricas_train <- bind_rows(metricas_mlm,metricas_semilog,
metricas_lmPropio1,metricas_lmPropio2)
metricas_train_horizontal <- pivot_wider(metricas_train,names_from=.metric,values_from=.estimate,id_cols=modelo)
metricas_train_horizontal$dataset <- "Train"eph_test <- eph_test %>% mutate(exp2=experiencia_potencial**2,
zona = case_when(region == "Capital Federal" ~ "AMBA",
region == "Gran Buenos Aires" ~ "AMBA",
region == "Noreste" ~ "Norte",
region == "Noroeste" ~ "Norte",
region == "Patagonia" ~ "Sur",
region == "Cuyo" ~ "Centro",
region == "Pampeana" ~ "Centro"),
vida_trabajada = experiencia_potencial/edad
)
# Obtenemos las predicciones de ambos modelos
lista_predicciones_testing = map(.x = models, .f = broom::augment, newdata = eph_test)
metricas_mlm = lista_predicciones_testing$mlm %>%
metrics(truth=salario_horario, estimate=.fitted) %>%
mutate(.estimate=round(.estimate, 4))
metricas_mlm$modelo <- "mlm"
metricas_lmPropio1 = lista_predicciones_testing$lmPropio1%>%
metrics(truth=salario_horario, estimate=.fitted) %>%
mutate(.estimate=round(.estimate, 4))
metricas_lmPropio1$modelo <- "lmPropio1"
metricas_lmPropio2 = lista_predicciones_testing$lmPropio2 %>%
metrics(truth=salario_horario, estimate=.fitted) %>%
mutate(.estimate=round(.estimate, 4))
metricas_lmPropio2$modelo <- "lmPropio2"
metricas_semilog = lista_predicciones_testing$lmMincerEnriquecido %>%
mutate(fitted_antilog= exp(.fitted)) %>%
metrics(truth=salario_horario, estimate=fitted_antilog) %>%
mutate(.estimate=round(.estimate, 4))
metricas_semilog$modelo <- "lmMincerEnriquecido"
#Juntamos las métricas de cada modelo
metricas <- bind_rows(metricas_mlm,metricas_semilog,metricas_lmPropio1,metricas_lmPropio2)
metricas_horizontal <- pivot_wider(metricas,names_from=.metric,values_from=.estimate,id_cols=modelo)
metricas_horizontal$dataset <- "Test"
metricas <- rbind(metricas_train_horizontal,metricas_horizontal)
metricas## # A tibble: 8 × 5
## modelo rmse rsq mae dataset
## <chr> <dbl> <dbl> <dbl> <chr>
## 1 mlm 376. 0.163 243. Train
## 2 lmMincerEnriquecido 386. 0.179 232. Train
## 3 lmPropio1 368. 0.202 234. Train
## 4 lmPropio2 376. 0.164 243. Train
## 5 mlm 385. 0.158 249. Test
## 6 lmMincerEnriquecido 397. 0.172 240. Test
## 7 lmPropio1 375. 0.200 240. Test
## 8 lmPropio2 386. 0.154 250. TestCon estos resultados, podemos afirmar que el modelo de mincer enriquecido tiene el minimo error absoluto y el modelo que considera la segmentación por zonas tiene la mejór performance en error cuadrático.
En este punto vale destacar que en caso de existir valores atípicos, la performance medida con error absoluto se verá menos afectada que con el error cuadrático. Si tuviera que elegir un modelo para predecir, utilizaría el semielástico de Mincer enriquecido.
6) Modelo lineal robusto
Leer el archivo “eph_train_outliers_2022.csv”. Este último consiste en el dataset original de train con la incorporación de algunas observaciones adicionales que pueden incluir valores atípicos.
Realizar dos gráficos del salario horario, uno para el dataset de entrenamiento sin outliers y otro para el dataset con outliers que permitan observar claramente la diferencia entre ambos sets de datos.
En primera instancia cargamos los datos con outliers.
eph_train_outliers <- read.csv("eph_train_outliers_2022.csv") %>% select(-c(codigo_actividad,codusu,aglomerado,ano4,trimestre)) %>% mutate(sexo=as.factor(sexo))Luego combinamos de manera vertical las variables de interes (edad, educacion, experiencia potencial y salario horario) para poder graficarlas de manera adecuada y observar las caracteristicas de los outliers.
eph_train2 <- eph_train %>% select(-c(exp2,zona,vida_trabajada))
eph_train_outliers2 <- eph_train_outliers
eph_train2$type <- "WO_outliers"
eph_train_outliers2$type <- "outliers"
combined <- rbind(eph_train2,eph_train_outliers2)
rm(eph_train2,eph_train_outliers2)
combinedPivot <- combined %>% select(where(is.numeric),type) %>% pivot_longer(.,cols=-c("type"))Ahora generamos graficos de cajas y bigotes de las distintas variables y las abrimos segun pertenezcan al dataset limpio o aquel que contiene outliers.
#Graficos univariados para ver outliers
ggplot(data = combinedPivot, aes(y = value, color = type)) +
geom_boxplot(alpha = 0.75) + # agregamos transparencia a los puntos
labs(title = "Variables numéricas de interés") +
facet_wrap(~ name, scales = "free")
Podemos ver que los outliers aparecen en valores altos de salario_horario. Ahora veamos cómo se ubican respecto de las demás variables numéricas, mediante gráficos de dispersion.
#Graficos bivariados para ver cómo podrian afectar los outliers.
combinedPivot2 <- combined %>% select(where(is.numeric),type) %>% pivot_longer(.,cols=-c("type","salario_horario")) %>% select(c(salario_horario,name,value,type))
ggplot(data = combinedPivot2, aes(x=value,y = salario_horario, color = type)) +
geom_point(alpha=0.2) + # agregamos transparencia a los puntos
labs(title = "Salario horario en función de otras variables numéricas") +
facet_wrap(~ name, scales = "free")
Los outliers de salario horario se corresponden con valores de edad, educacion y experiencia potencial distribuidos en todo el ámbito de estas variables. Esto quiere decir que no necesariamente todos los llamados “outliers” tendrán gran capacidad de apalancamiento.
Sobre este nuevo conjunto de datos entrenar el modelo lineal multiple, el modelo de mincer y un modelo robusto (misma especificación que el modelo lineal multiple). Comparar exhaustivamente los coeficientes estimados y su significatividad entre el modelo lineal multiple y el modelo robusto. Comparar la performance (RMSE y MAE) de los tres modelos entrenados en este punto en el dataset de entrenamiento (con outliers) y de evaluación ¿Qué puede concluir al respecto?
#Entrenamos modelos con outliers
mlm_outliers <- eph_train_outliers %>% mutate(exp2=experiencia_potencial**2) %>% lm(salario_horario~educacion+experiencia_potencial+exp2+sexo+sexo*educacion,.)
summary(mlm_outliers)##
## Call:
## lm(formula = salario_horario ~ educacion + experiencia_potencial +
## exp2 + sexo + sexo * educacion, data = .)
##
## Residuals:
## Min 1Q Median 3Q Max
## -950.5 -223.5 -76.0 105.2 28627.2
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -324.05138 36.09441 -8.978 < 2e-16 ***
## educacion 49.24487 2.24080 21.976 < 2e-16 ***
## experiencia_potencial 10.18976 1.31632 7.741 1.07e-14 ***
## exp2 -0.07105 0.02502 -2.840 0.00452 **
## sexoVaron 52.01762 39.97978 1.301 0.19325
## educacion:sexoVaron 0.26172 2.90051 0.090 0.92810
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 558.4 on 11636 degrees of freedom
## Multiple R-squared: 0.0881, Adjusted R-squared: 0.08771
## F-statistic: 224.8 on 5 and 11636 DF, p-value: < 2.2e-16lmMincerEnriquecido_outliers <- eph_train_outliers %>% mutate(exp2=experiencia_potencial**2) %>% lm(log(salario_horario) ~ educacion+experiencia_potencial+exp2 + sexo +
sexo*educacion,.)
summary(lmMincerEnriquecido_outliers)##
## Call:
## lm(formula = log(salario_horario) ~ educacion + experiencia_potencial +
## exp2 + sexo + sexo * educacion, data = .)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.5973 -0.3665 0.0354 0.3919 4.2231
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.4405110 0.0414047 107.247 < 2e-16 ***
## educacion 0.0899883 0.0025705 35.008 < 2e-16 ***
## experiencia_potencial 0.0232941 0.0015100 15.427 < 2e-16 ***
## exp2 -0.0002656 0.0000287 -9.257 < 2e-16 ***
## sexoVaron 0.2387073 0.0458616 5.205 1.97e-07 ***
## educacion:sexoVaron -0.0106225 0.0033272 -3.193 0.00141 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6406 on 11636 degrees of freedom
## Multiple R-squared: 0.1807, Adjusted R-squared: 0.1803
## F-statistic: 513.2 on 5 and 11636 DF, p-value: < 2.2e-16#Usamos el modelo robusto lmRob
mlm_robusto_outliers <- eph_train_outliers %>% mutate(exp2=experiencia_potencial**2) %>% lmRob(salario_horario~educacion+experiencia_potencial+exp2 +sexo+sexo*educacion,.)
summary(mlm_robusto_outliers)##
## Call:
## lmRob(formula = salario_horario ~ educacion + experiencia_potencial +
## exp2 + sexo + sexo * educacion, data = .)
##
## Residuals:
## Min 1Q Median 3Q Max
## -566.5491 -140.7266 -0.9109 184.3904 28715.7651
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -150.04098 16.12079 -9.307 < 2e-16 ***
## educacion 34.31684 1.01121 33.936 < 2e-16 ***
## experiencia_potencial 8.47156 0.58458 14.492 < 2e-16 ***
## exp2 -0.09689 0.01117 -8.671 < 2e-16 ***
## sexoVaron 102.02343 17.87848 5.706 1.18e-08 ***
## educacion:sexoVaron -5.62929 1.31469 -4.282 1.87e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 240.4 on 11636 degrees of freedom
## Multiple R-Squared: 0.1116
##
## Test for Bias:
## statistic p-value
## M-estimate -297.1 1
## LS-estimate 3916.7 0Comparación de los valores de los parámetros estimados
La ordenada al origen para las mujeres es entre 2 y 3 veces más grande en módulo en el modelo lineal múltiple que en el robusto. El coeficiente que acompaña a educacion (en el caso de mujeres) es mayor (entre 1.5 y 2 veces) que en el modelo robusto. En los coeficientes que acompañan experiencia potencial y su cuadrado, las estimaciones con ambas estrategias no son distinguibles (los errores estándar se solapan). En el termino independiente que atribuye la diferencia entre Varones y Mujeres, el método robusto duplica en módulo al coeficiente del modelo comun. Por último, la diferencia en la influencia de la educacion en varones sobre el salario horario esperado, por sobre las mujeres, tiene signo opuesto y es aproximadamente 20 veces mayor en módulo en el modelo robusto que en el lineal comun.
Si hacemos un poco de zoom en los valores donde el ajuste es mejor y graficamos los modelos lineales usando educacion como covariable
Ploteamos los dos modelos Nivel-Nivel en funcion de educacion
#Buscamos coeficientes para educacion y graficamos las rectas
intercept_mlm_mujer=mlm_outliers$coefficients[1]
slope_mlm_mujer=mlm_outliers$coefficients[2]
intercept_mlm_varon=mlm_outliers$coefficients[1]+mlm_outliers$coefficients[5]
slope_mlm_varon=mlm_outliers$coefficients[2]+mlm_outliers$coefficients[6]
tidy(mlm_robusto_outliers)## # A tibble: 6 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) -150. 16.1 -9.31 1.55e- 20
## 2 educacion 34.3 1.01 33.9 8.41e-241
## 3 experiencia_potencial 8.47 0.585 14.5 3.52e- 47
## 4 exp2 -0.0969 0.0112 -8.67 4.85e- 18
## 5 sexoVaron 102. 17.9 5.71 1.18e- 8
## 6 educacion:sexoVaron -5.63 1.31 -4.28 1.87e- 5intercept_rmlm_mujer=mlm_robusto_outliers$coefficients[1]
slope_rmlm_mujer=mlm_robusto_outliers$coefficients[2]
intercept_rmlm_varon=mlm_robusto_outliers$coefficients[1]+mlm_robusto_outliers$coefficients[5]
slope_rmlm_varon=mlm_robusto_outliers$coefficients[2]+mlm_robusto_outliers$coefficients[6]
eph_train_outliers %>% ggplot(., aes(x = educacion, y = salario_horario)) +
geom_abline(intercept = intercept_mlm_mujer, slope = slope_mlm_mujer, color="red") +
geom_abline(intercept = intercept_mlm_varon, slope = slope_mlm_varon, color="green") +
geom_abline(intercept = intercept_rmlm_mujer, slope = slope_rmlm_mujer, color="magenta") +
geom_abline(intercept = intercept_rmlm_varon, slope = slope_rmlm_varon, color="forestgreen") +
geom_point() +
ylim(0, 2000)+
xlim(0, 30)+
labs(title="Modelo Lineal Simple", x="Educacion(años)", y="Salario Horario (ars/h)") ## Warning: Removed 158 rows containing missing values (`geom_point()`).
Si comparamos los valores de los parámetros estimados en la regresion lineal múltiple sobre los datasets con y sin outliers veremos que no se modifican significativamente los coeficientes relacionados con educacion y experiencia potencial, mientras que los relacionados con el sexo y la interacción entre sexo y educación varian de manera apreciable.
Respecto del modelo robusto(curvas fucsia y verde oscuro), el modelo múltiple comun tiene mayor pendiente y menor ordenada al origen (curvas roja y verde clara en la figura anterior), dando a entender que la asimetria en la distribución de salarios para una dada educacion, experiencia laboral y sexo, tironean de este parámetro haciendolo variar (en términos absolutos) entre 1.5 y 3 veces más. Consecuentemente, esto hace que la ordenada del modelo lineal comun sea un valor negativo de mayor módulo que el obtenido en el modelo robusto. El sexo como variable categorica también modifica su influencia sobre el salario horario medio en un factor de 2, pasando de ser de aproximadamente 52 ars/h en el modelo comun y 100 en el robusto, a igualdad del resto de las covariables.
Por otro lado, resulta interesante ver qué ocurre con el término de interacción: en el caso del modelo lineal común, el efecto de interacción es muy pequeño y los modelos resultantes para Mujeres y Varones resultan prácticamente paralelos en todo el rango de covariables estudiado. Sin embargo, en el caso robusto, puede verse que este coeficiente (educacion:sexoVaron) tiene un valor negativo. Esto quiere decir que un incremento de la variable educación repercute sobre los Varones de manera negativa sobre la estimacion informada para el sexo Mujer (o sea, al ser este coeficiente estimado como -5.6 ars/( h * año educacion), estamos afirmando que el coeficiente de educacion para la clase de referencia (33.2 ars/( h * año educacion)) disminuye en esa cantidad, dando un coeficiente efectivo para educacion en el grupo de los Varones de 27.6 ars/(h * año educacion).
Significatividad de estos coeficientes
Para el modelo robusto, todos los coeficientes obtenidos resultan significativos a un nivel de confianza de hasta 0.9998. Por otro lado, en el modelo lineal multiple original, los términos que involucran la distincion entre sexos (término independiente y modificacion al término que refiere a educación) tienen p valores muy altos (mayores a 0.2), lo cual habla sobre la poca evidencia que puede ponderar dicho modelo para afirmar que estos parámetros son distintos de cero.
Métricas de performance para los tres modelos
Estimamos RMSE y MAE para el conjunto de train (con outliers) y test.
En train
models <- list(mlm_outliers=mlm_outliers,lmMincerEnriquecido_outliers = lmMincerEnriquecido_outliers, mlm_robusto_outliers=mlm_robusto_outliers)
eph_train_outliers_complete <- eph_train_outliers %>% mutate(exp2=experiencia_potencial**2,
zona = case_when(region == "Capital Federal" ~ "AMBA",
region == "Gran Buenos Aires" ~ "AMBA",
region == "Noreste" ~ "Norte",
region == "Noroeste" ~ "Norte",
region == "Patagonia" ~ "Sur",
region == "Cuyo" ~ "Centro",
region == "Pampeana" ~ "Centro"),
vida_trabajada = experiencia_potencial/edad,
sexo=as.factor(sexo))
lista_predicciones_train = map(.x = models, .f = broom::augment,newdata=eph_train_outliers_complete)
metricas_mlm_outliers = lista_predicciones_train$mlm_outliers %>% as_tibble() %>%
metrics(truth=salario_horario, estimate=.fitted) %>%
mutate(.estimate=round(.estimate, 4))
metricas_mlm_outliers$modelo <- "mlm_outliers"
metricas_semilog_outliers = lista_predicciones_train$lmMincerEnriquecido_outliers %>% as_tibble() %>% mutate(fitted_antilog= exp(.fitted)) %>% metrics(truth=salario_horario, estimate=fitted_antilog) %>%
mutate(.estimate=round(.estimate, 4))
metricas_semilog_outliers$modelo <- "lmMincerEnriquecido_outliers"
metricas_mlm_rob_outliers = lista_predicciones_train$mlm_robusto_outliers %>%
metrics(truth=salario_horario, estimate=.fitted) %>%
mutate(.estimate=round(.estimate, 4))
metricas_mlm_rob_outliers$modelo <- "mlm_robusto_outliers"
metricas <- bind_rows(metricas_mlm_outliers,metricas_semilog_outliers,metricas_mlm_rob_outliers)
metricas_horizontal <- pivot_wider(metricas,names_from=.metric,values_from=.estimate,id_cols=modelo)
metricas_horizontal$datos <- "Train"En test
eph_test <- eph_test %>% mutate(sexo=as.factor(sexo))
lista_predicciones_test = map(.x = models, .f = broom::augment,newdata=eph_test)
metricas_mlm_outliers2 = lista_predicciones_test$mlm_outliers %>% as_tibble() %>%
metrics(truth=salario_horario, estimate=.fitted) %>%
mutate(.estimate=round(.estimate, 4))
metricas_mlm_outliers2$modelo <- "mlm_outliers"
metricas_semilog_outliers2 = lista_predicciones_test$lmMincerEnriquecido_outliers %>% as_tibble() %>% mutate(fitted_antilog= exp(.fitted)) %>% metrics(truth=salario_horario, estimate=fitted_antilog) %>%
mutate(.estimate=round(.estimate, 4))
metricas_semilog_outliers2$modelo <- "lmMincerEnriquecido_outliers"
metricas_mlm_rob_outliers2 = lista_predicciones_test$mlm_robusto_outliers %>%
metrics(truth=salario_horario, estimate=.fitted) %>%
mutate(.estimate=round(.estimate, 4))
metricas_mlm_rob_outliers2$modelo <- "mlm_robusto_outliers"
metricas_test <- bind_rows(metricas_mlm_outliers2,metricas_semilog_outliers2,metricas_mlm_rob_outliers2)
metricas_test_horizontal <- pivot_wider(metricas_test,names_from=.metric,values_from=.estimate,id_cols=modelo)
metricas_test_horizontal$datos <- "Test"
metricas <- rbind(metricas_horizontal,metricas_test_horizontal)
metricas## # A tibble: 6 × 5
## modelo rmse rsq mae datos
## <chr> <dbl> <dbl> <dbl> <chr>
## 1 mlm_outliers 558. 0.0881 260. Train
## 2 lmMincerEnriquecido_outliers 567. 0.0953 245. Train
## 3 mlm_robusto_outliers 570. 0.0859 248. Train
## 4 mlm_outliers 385. 0.157 252. Test
## 5 lmMincerEnriquecido_outliers 396. 0.172 240. Test
## 6 mlm_robusto_outliers 399. 0.158 241. TestLas métricas de performance RMSE y MAE tanto en los datasets de entrenamiento (con outliers) y de testeo reconocen una mejor performance en RMSE al modelo lineal múltiple y al modelo de mincer enriquecido. Es notable el incremento (casi al doble) de la raiz del error cuadrático en datasets de entrenamiento respecto de la validación. Esto puede deberse a la gran dispersión que tienen los datos en entrenamiento, y que el aporte que tienen valores dispersos en la métrica cuadrática es mayor.
Por otro lado, el error absoluto medio es semejante al evaluar en entrenamiento y validación, al poderar linealmente las distancias al valor medio. Los extremos afectan de una forma más moderada esta métrica.
Si nos guiamos por la métrica MAE, predecir utilizando el modelo de Mincer enriquecido es el que otorgará en promedio mejores predicciones. Surge la pregunta: Pareceria ser que el ajuste de un modelo robusto no trajo mejoras significativas para predecir nuevos valores ¿por qué?
Una hipótesis sería que los outliers eliminados por los docentes sean realmente datos atípicos pero que no tienen un potencial de apalancamiento grande. Para esto,observamos el gráfico de residuos versus apalancamiento potencial.
#Obtenemos los graficos para realizar el diagnóstico
plot(mlm_outliers)
Podemos observar que todos los puntos considerados se encuentran en la región donde los datos no representan un riesgo a la hora de afectar el ajuste y la capacidad de predicción del modelo (todos los puntos se encuentran dentro de la región delimitada por la distancia de Cook 0.5). En caso de aparecer puntos con distancia de Cook mayor a 0.5, me tomaria el trabajo de evaluar nuevamente la performance relativa de estos modelos, dado que en esa instancia el modelo robusto podria traer ventajas en la predicción de nuevos valores.
Bibliografia consultada
Apuntes de la materia
Apunte de Regresión Lineal - María Eugenia Szretter Noste
Documentación online de las librerias de R.